Quelques explorations cartographiques avec les données OSM
Timothée Giraud
17 novembre 2017
1 Téléchargement et préparation des données
1.1 Téléchargement des données
library(osmdata)
library(sf)
library(units)
# define a bounding box
q0 <- opq(bbox = c( 2.2247, 48.8188, 2.4611, 48.9019))
# extract Paris boundaries
q1 <- add_osm_feature(opq = q0, key = 'name', value = "Paris", value_exact = TRUE)
res1 <- osmdata_sf(q1)
paris <- st_geometry(res1$osm_multipolygons[1,])
# extract La Seine
q2 <- add_osm_feature(opq = q0, key = 'name', value = "La Seine", value_exact = TRUE)
res2 <- osmdata_sf(q2)
seine1 <- st_geometry(res2$osm_multilines)
q2b <- add_osm_feature(opq = q0, key = 'name', value = "La Seine - Bras de la Monnaie", value_exact = FALSE)
res2b <- osmdata_sf(q2b)
seine2 <- st_geometry(res2b$osm_lines)
# extract Parcs and Cimetierres
q3 <- add_osm_feature(opq = q0, key = 'leisure', value = "park", value_exact = TRUE)
res3 <- osmdata_sf(q3)
parc1 <- st_geometry(res3$osm_polygons)
parc2 <- st_geometry(res3$osm_multipolygons)
q4 <- add_osm_feature(opq = q0, key = 'landuse', value = "cemetery", value_exact = TRUE)
res4 <- osmdata_sf(q4)
parc3 <- st_geometry(res4$osm_polygons)
# extract Quartiers
q5 <- add_osm_feature(opq = q0, key = 'admin_level', value = "10", value_exact = TRUE)
res5 <- osmdata_sf(q5)
quartier <- res5$osm_multipolygons
# extract Bars, Pubs & Cafe
q6 <- add_osm_feature(opq = q0, key = 'amenity', value = "bar", value_exact = TRUE)
res6 <- osmdata_sf(q6)
bar <- res6$osm_points
q7 <- add_osm_feature(opq = q0, key = 'amenity', value = "cafe", value_exact = TRUE)
res7 <- osmdata_sf(q7)
cafe <- res7$osm_points
q8 <- add_osm_feature(opq = q0, key = 'amenity', value = "pub", value_exact = TRUE)
res8 <- osmdata_sf(q8)
pub <- res8$osm_points
# extract Metro Stations
q9 <- add_osm_feature(opq = q0, key = 'station', value = "subway", value_exact = TRUE)
res9 <- osmdata_sf(q9)
metro <- res9$osm_points
# extract Restaurant
q10 <- add_osm_feature(opq = q0, key = 'amenity', value = "restaurant", value_exact = TRUE)
res10 <- osmdata_sf(q10)
resto <- res10$osm_points1.2 Préparation des données
# use Lambert 93 projection
parc1 <- st_transform(parc1, 2154)
parc2 <- st_transform(parc2, 2154)
parc3 <- st_transform(parc3, 2154)
paris <- st_transform(paris, 2154)
seine1 <- st_transform(seine1, 2154)
seine2 <- st_transform(seine2, 2154)
quartier <- st_transform(quartier, 2154)
bar <- st_transform(bar, 2154)
pub <- st_transform(pub, 2154)
cafe <- st_transform(cafe, 2154)
metro <- st_transform(metro, 2154)
resto <- st_transform(resto, 2154)
# make layers pretty
parc <- do.call(c, list(parc1, parc2, parc3))
parc <- st_union(x = st_buffer(parc,0), by_feature = F)
parc <- st_cast(parc, "POLYGON")
parc <- parc[st_area(parc)>=set_units(10000, "m^2")]
parc <- st_intersection(x = parc, y = paris)
seine <- st_intersection(x = seine1, y = paris)
seine <- c(st_cast(seine[1])[2:5], seine[2])
seine <-c(seine, seine2)
bar <- bar[!is.na(bar$name),]
pub <- pub[!is.na(pub$name),]
cafe <- cafe[!is.na(cafe$name),]
bars <- c(st_geometry(bar), st_geometry(pub))
# bars <- c(st_geometry(bar), st_geometry(cafe), st_geometry(pub))
bars <- st_intersection(x = bars, y = paris)
quartier <- quartier[substr(quartier$ref.INSEE,1,2)==75,]
metro <- st_intersection(x = metro, y = paris)
metro <- metro[metro$station%in%"subway", ]
metro <- metro[metro$railway%in%"station",]
metro <- metro[metro$STIF.zone%in%1,]
resto <- resto[resto$amenity%in%"restaurant",]
resto <- resto[!is.na(resto$name),]
resto <- resto[,"cuisine"]
resto <- st_intersection(x = resto, y = paris)
resto$cuisine <- as.character(resto$cuisine)
save(list= c("paris", "quartier", "seine", "parc", "bars", "metro", "resto"),
file = "../data/paname.RData", compress = "xz")2 Analyse et Cartographie des bars à Paris
2.1 Affichage des données
library(spatstat)
library(sf)
library(cartography)
library(maptools)
library(raster)
library(units)
load("../data/paname.RData")
# 1st map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
plot(bars, add=T, col = "#0000ff", pch = 20, cex = 0.2)
layoutLayer(title = "Les bars à Paris", scale = 1,
north = T,
tabtitle = TRUE,
sources = "Map data © OpenStreetMap contributors, under CC BY SA.",
author = "cartography 2.0.2") 
2.2 Nombre de bars par quartiers
# compute some data
quartier$nbar <- lengths(st_covers(quartier, bars))
quartier$dbar <- quartier$nbar / set_units(st_area(quartier), "km^2")
# 2nd map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
propSymbolsLayer(x = quartier, var = "nbar", inches = 0.2, col = "red",
legend.pos = "topright",
legend.title.txt = "Nombre de bars\npar quartier")
layoutLayer(title = "Les bars à Paris", scale = 1,
tabtitle = TRUE,
sources = "Map data © OpenStreetMap contributors, under CC BY SA.",
author = "cartography 2.0.2") 
# 3rd map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
choroLayer(quartier, var = "dbar", border = NA,
method = "quantile", nclass = 6,
col = carto.pal("wine.pal", 6),
legend.pos = "topright",
legend.title.txt = "Densité de bars\npar km2", add = TRUE)
plot(parc, col = "#E2CCB5", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
layoutLayer(title = "Les bars à Paris", scale = 1,
tabtitle = TRUE,
sources = "Map data © OpenStreetMap contributors, under CC BY SA.",
author = "cartography 2.0.2") 
# 2nd map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
propSymbolsChoroLayer(x = quartier, var = "nbar", var2 = "dbar",
inches = 0.2,
method = "quantile", nclass = 6,
col = carto.pal("wine.pal", 6),
legend.var2.pos = "topright",
legend.var2.title.txt = "Densité de bars\npar km2",
legend.var.title.txt = "Nombre de bars\npar quartier")
layoutLayer(title = "Les bars à Paris", scale = 1,
tabtitle = TRUE,
sources = "Map data © OpenStreetMap contributors, under CC BY SA.",
author = "cartography 2.0.2") 
2.3 Analyse par station de métro
library(sf)
# distance entre toutes les stations
x <- st_distance(metro, metro)
diag(x) <- 10000
# distance à la plus proche station
dmin <- apply(x, 2, min)
metro$nbar <- lengths(st_covers(st_buffer(x = metro, dist = mean(dmin)), bars))
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
propSymbolsChoroLayer(x = metro, var = "nbar", var2 = "nbar",
col = carto.pal("wine.pal",5), inches = 0.1,
method = "fisher-jenks", nclass = 5, add=T,
legend.var.pos = "topleft",
legend.var.title.txt = "Nombre de bars à proximité\nde la station de métro")
layoutLayer(title = "Bars à Paris", scale = 1,
tabtitle = TRUE,
sources = "Map data © OpenStreetMap contributors, under CC BY SA.",
author = "cartography 2.0.2") 
2.4 Carte interactive
library(leaflet)
## Initialisation
metrounp <- st_transform(metro, 4326)
metrounp <- metrounp[order(metrounp$nbar, decreasing = T),]
metrounp$size <- sqrt(metrounp$nbar*15 / pi)
metrounp$label <- paste0("<b>", metrounp$name, "</b> <br>",
metrounp$nbar," bars à proximité.")
m <- leaflet(padding = 0)
m <- addTiles(m)
m <- addCircleMarkers(map = m,
lng = st_coordinates(metrounp)[,1],
lat = st_coordinates(metrounp)[,2],
radius = metrounp$size, weight = 0.25,
stroke = T, opacity = 100,
fill = T, fillColor = "#920000",
fillOpacity = 100,
popup = metrounp$label,
color = "white")
m2.5 Analyse lissée de la répartition des bars
bb <- as(bars, "Spatial")
bbowin <- as.owin(as(paris, "Spatial"))
pts <- coordinates(bb)
p <- ppp(pts[,1], pts[,2], window=bbowin)
ds <- density.ppp(p, sigma = 200, eps = c(20,20))
rasdens <- raster(ds) * 1000 * 1000
rasdens <- rasdens+1
par(mar = c(0,0,1.2,0))
bks <- getBreaks(values(rasdens), nclass = 12, method = "equal")
cols <- colorRampPalette(c("black","#940000", "white"))(length(bks)-1)
plot(paris, col = NA, border = NA, main="", bg = "#FBEDDA")
plot(rasdens, breaks= bks, col=cols, add = T,legend=F)
legendChoro(pos = "topright",cex = 0.7,
title.txt = "Densité de Bars\nkde, sigma=200m\n(bars par km2)",
breaks = bks-1, nodata = FALSE,
col = cols)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(parc, col = "#CDEBB235", border = NA, add=T)
plot(quartier$geometry, add=T, border = "white",lty = 2, lwd = 0.1)
plot(bars, add=T, col = "#0000ff90", pch = 20, cex = 0.1)
plot(paris, col = NA, add=T)
layoutLayer(title = "Les bars à Paris", scale = 1,
tabtitle = TRUE,
sources = "Map data © OpenStreetMap contributors, under CC BY SA.",
author = "cartography 2.0.2") 
3 Analyse et cartographie des restaurants à Paris
3.1 Les données
# 1st map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
plot(resto, add=T, col = "blue", pch = 20, cex = 0.1)
layoutLayer(title = "Tous les restaurants", scale = 1,
north = T,
tabtitle = TRUE,
sources = "Map data © OpenStreetMap contributors, under CC BY SA.",
author = "cartography 2.0.2") 
3.2 Cartes de densités
densityMap <- function(x = resto, cuisine = NULL, title, sigma ){
opar <- par(mar = c(0,0,0,0))
if(!is.null(cuisine)){
x <- x[x$cuisine %in% cuisine, ]
}
bb <- as(x, "Spatial")
bbowin <- as.owin(as(paris, "Spatial"))
pts <- coordinates(bb)
p <- ppp(pts[,1], pts[,2], window=bbowin)
ds <- density.ppp(p, sigma = sigma, eps = c(20,20))
rasdens <- raster(ds) * 1000 * 1000
rasdens <- rasdens+1
xx <- getBreaks(values(rasdens), nclass = 12, method = "equal")
cols <- colorRampPalette(c("black", "#940000", "white"))(length(xx)-1)
plot(paris, col = NA, border = NA, main="")
image(rasdens, breaks= xx, col=cols, add = T,legend=F)
legendChoro(pos = "topright",cex = 0.7,title.cex = 1,
title.txt = paste0(title, "\nn=",nrow(pts)),
breaks = xx-1, nodata = FALSE,values.rnd = 0,
col = cols)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(parc, col = "#CDEBB235", border = NA, add=T)
plot(quartier$geometry, add=T, border = "white",lty = 2, lwd = 0.1)
plot(st_geometry(x), add=T, col = "blue", pch = 20, cex = 0.1)
barscale(size = 1)
mtext(text = "cartography 2.0.2\nMap data © OpenStreetMap contributors, under CC BY SA.",
side = 1, line = -1, adj = c(0.01,0.01), cex = 0.6, font = 3)
north(pos = c(661171.8,6858051))
par(opar)
}
densityMap(x = resto, title = "Tous les restaurants", sigma = 500)
densityMap(x = resto, cuisine = c('japanese','sushi'), title = "Cuisine japonaise", sigma = 500)
densityMap(x = resto, cuisine = c('pizza','italian'), title = "Cuisine italienne", sigma = 500)
densityMap(x = resto, cuisine = c('chinese'), title = "Cuisine chinoise", sigma = 500)
densityMap(x = resto, cuisine = c('korean'), title = "Cuisine coréenne", sigma = 500)
densityMap(x = resto, cuisine = c('indian'), title = "Cuisine Indienne", sigma = 500)
densityMap(x = resto, cuisine = c('asian', "thai", "vietnamese"), title = "Cuisine asiatique", sigma = 500)
4 Pour la reproductibilité…
Voici les informations sur ma configuration :
sessionInfo()## R version 3.4.2 (2017-09-28)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## locale:
## [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
## [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
## [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] leaflet_1.1.0 units_0.4-6 raster_2.6-7
## [4] maptools_0.9-2 cartography_2.0.2 sp_1.2-5
## [7] sf_0.5-5 spatstat_1.53-2 rpart_4.1-11
## [10] nlme_3.1-131 spatstat.data_1.2-0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.13 compiler_3.4.2 class_7.3-14
## [4] tools_3.4.2 digest_0.6.12 goftest_1.1-1
## [7] jsonlite_1.5 evaluate_0.10.1 lattice_0.20-35
## [10] mgcv_1.8-22 Matrix_1.2-11 shiny_1.0.5
## [13] DBI_0.7 crosstalk_1.0.0 yaml_2.1.14
## [16] e1071_1.6-8 stringr_1.2.0 knitr_1.17
## [19] htmlwidgets_0.9 rgeos_0.3-26 spatstat.utils_1.8-0
## [22] classInt_0.1-24 rprojroot_1.2 grid_3.4.2
## [25] R6_2.2.2 foreign_0.8-69 rmarkdown_1.8
## [28] polyclip_1.6-1 deldir_0.1-14 udunits2_0.13
## [31] magrittr_1.5 tensor_1.5 backports_1.1.1
## [34] htmltools_0.3.6 abind_1.4-5 xtable_1.8-2
## [37] mime_0.5 httpuv_1.3.5 stringi_1.1.6